Load required packages

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.1     ✔ purrr   1.0.1
## ✔ tibble  3.1.8     ✔ dplyr   1.1.0
## ✔ tidyr   1.3.0     ✔ stringr 1.5.0
## ✔ readr   2.1.4     ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(stringi)
library(doParallel)
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## 
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## 
## Loading required package: iterators
## Loading required package: parallel
library(pomp)
## 
## Welcome to pomp!
## 
## As of version 4.6, no user-visible pomp function has a name that
## includes a dot ('.'). Function names have been changed to replace the
## dot with an underscore ('_'). For more information, see the pomp blog:
## https://kingaa.github.io/pomp/blog.html.
## 
## 
## Attaching package: 'pomp'
## 
## The following object is masked from 'package:purrr':
## 
##     map
library(panelPomp)
library(ggpubr)

Set working directory

setwd("~/Documents/GitHub/bdd/nw11_hier/final_files/")

Colour-blind friendly colour palette

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

Read in PNAS data

## # A tibble: 6 × 11
##     day       Pd     RBC   Ter119   CD71 mouseid paba  box   mouse  Eryth  Retic
##   <int>    <dbl>   <dbl>    <dbl>  <dbl> <chr>   <chr> <chr> <chr>  <dbl>  <dbl>
## 1     0      NA  8360000 5105135. 1.37e5 01-01   0.05  01    01    8.14e6 2.25e5
## 2     1    7934. 8290000 2745516. 1.94e5 01-01   0.05  01    01    7.70e6 5.86e5
## 3     2   19489. 7560000 2551937. 1.51e5 01-01   0.05  01    01    7.11e6 4.46e5
## 4     3  228842. 7820000 6400210. 4.98e5 01-01   0.05  01    01    7.21e6 6.08e5
## 5     4 1534425  7520000 3471975. 4.16e5 01-01   0.05  01    01    6.62e6 9.01e5
## 6     5 4528560  7600000 4748225. 5.50e5 01-01   0.05  01    01    6.72e6 8.80e5

Obtain estimates for beta and dose

##   box mouse mouseid     Beta       dose sigmaPd sigmaRBC sigmaRetic sigmaW
## 1  01    01   01-01 6.511778 793.298642       2      0.1        0.3      1
## 2  01    02   01-02 6.511778   3.444740       2      0.1        0.3      1
## 3  01    03   01-03 6.511778 420.158396       2      0.1        0.3      1
## 4  02    01   02-01 6.090502 607.635623       2      0.1        0.3      1
## 5  02    02   02-02 6.090502 260.448444       2      0.1        0.3      1
## 6  02    03   02-03 6.090502   1.348118       2      0.1        0.3      1
##   sigmaN sigmaR   E_0   R_0   W_0   N_0
## 1    0.5    0.5 8e+06 3e+05 88000 8e+05
## 2    0.5    0.5 8e+06 3e+05 88000 8e+05
## 3    0.5    0.5 8e+06 3e+05 88000 8e+05
## 4    0.5    0.5 8e+06 3e+05 88000 8e+05
## 5    0.5    0.5 8e+06 3e+05 88000 8e+05
## 6    0.5    0.5 8e+06 3e+05 88000 8e+05

Create object pos with pomp object for each mouse

foreach (m = iter(theta,"row"),.inorder=TRUE,.combine=c) %dopar% {
  
  flow %>%
    filter(mouseid==m$mouseid) %>%
    select(day,Pd,RBC,Retic) %>%
    mutate(Retic=if_else(day %in% c(0,14),NA_real_,Retic)) %>%
    pomp(
      params=select(m,-mouseid,-box,-mouse) %>% unlist(),
      times="day",
      t0=0,
      rmeasure=Csnippet("
  Retic = rlnorm(log(1+R),sigmaRetic)-1;
  RBC = rlnorm(log(1+E+R),sigmaRBC)-1;
  Pd = rlnorm(log(1+K),sigmaPd)-1;"),
      dmeasure=Csnippet("
  double l1, l2, l3;
  l1 = (R_FINITE(Retic)) ? dlnorm(1+Retic,log(1+R),sigmaRetic,1) : 0;
  l2 = (R_FINITE(RBC)) ? dlnorm(1+RBC,log(1+E+R),sigmaRBC,1) : 0;
  l3 = (R_FINITE(Pd) && Pd>0) ? dlnorm(1+Pd,log(1+K),sigmaPd,1) : 0;
  lik = (give_log) ? l1+l2+l3 : exp(l1+l2+l3);"),
      rprocess=discrete_time(
        step.fun=Csnippet("
  double Mold = M;
  M = Beta*K*exp(-(W+N)/(R+E));
  E = (R+E)*exp(-(Mold+N)/(R+E));
  N = rlnorm(log(N),sigmaN);
  W = rlnorm(log(W),sigmaW);
  R = rlnorm(log(R),sigmaR);
  K = (R+E>0) ? (R+E)*(1-exp(-M/(R+E))): 0;
  "),
        delta.t=1
      ),
      partrans=parameter_trans(
        log=c("sigmaW","sigmaR","sigmaN",
              "sigmaPd","sigmaRBC","sigmaRetic",
              "N_0","W_0","E_0","R_0")
      ),
      rinit=Csnippet("
  E = E_0;
  R = R_0;
  N = N_0;
  W = W_0;
  M = 0;
  K = dose;"),
      statenames=c("E","R","W","N","M","K"),
      paramnames=c(
        "Beta","dose",
        "sigmaPd","sigmaRBC","sigmaRetic",
        "sigmaW","sigmaN","sigmaR",
        "E_0","R_0","W_0","N_0"
      )
    )
} %>% 
  set_names(theta$mouseid) -> pos

Obtain MLE for sigmas, initial values, betas and dose

##      loglik   sigmaW   sigmaN   sigmaPd   sigmaR   sigmaRBC sigmaRetic     E_0
## 1 -11217.79 1.043139 0.624558 0.5966612 0.392552 0.07942614  0.1514119 7722711
##        N_0      R_0      W_0
## 1 410608.7 617657.5 176481.2
## # A tibble: 15 × 3
##    mouseid  Beta    dose
##    <chr>   <dbl>   <dbl>
##  1 01-01    6.51  793.  
##  2 01-02    6.51    3.44
##  3 01-03    6.51  420.  
##  4 02-01    6.09  608.  
##  5 02-02    6.09  260.  
##  6 02-03    6.09    1.35
##  7 03-01    5.76  591.  
##  8 03-02    5.76  382.  
##  9 03-03    5.76  317.  
## 10 04-01    4.08 1476.  
## 11 04-02    4.08  716.  
## 12 04-03    4.08 1139.  
## 13 05-01    0       0   
## 14 05-02    0       0   
## 15 05-03    0       0

Load in PNAS trajectories

sm1name <- "m5sm1.rds"
sm1 <- readRDS(sm1name)

Create dataframe with weighted trajectories (grouped by mouseid first, then by box)

sm1 |>
  as_tibble() |>
  filter(mouseid!="01-02",mouseid!="02-03") |> #remove underdosed mice
  pivot_wider(names_from=variable,values_from=value) |>
  left_join(bdf,by=c("mouseid")) |>
  separate_wider_delim(cols="mouseid",delim="-",names=c("box","mouse"),cols_remove=FALSE) |>
  group_by(mouseid) |>
  mutate(
    SM=exp(-M/(R+E)),
    SN=exp(-N/(R+E)),
    SW=exp(-W/(R+E)),
    Qps=(1-SM)*SW*SN,
    Qpn=N/(N+W)*(1-SM)*(1-SW*SN),
    Qpw=W/(N+W)*(1-SM)*(1-SW*SN),
    Qun=SM*(1-SN),
    Qus=SM*SN,
    lambda_r=Beta*R/M*Qps,
    lambda_e=Beta*E/M*Qps,
    lambda_n=Beta*(R+E)/M*Qpn,
    lambda_w=Beta*(R+E)/M*Qpw,
    lambda_u=Beta-lambda_r-lambda_e-lambda_n-lambda_w,
    lambda_u_w_ratio=lambda_u/lambda_w,
    lambda_u_n_ratio=lambda_u/lambda_n,
    lambda_u_wn_ratio=lambda_u/(lambda_w+lambda_n),
    loss=(R+E)*(1-Qus),
    rbc=E+R,
    varN=N/rbc,
    perRetic=R/rbc,
    lik=exp(loglik-max(loglik))
  ) |>
  ungroup() |>
  select(-loglik,-Beta,-dose) |>
  gather(variable,value,-rep,-time,-mouse,-box,-lik,-mouseid) |>
  filter(is.finite(value)) %>%
  group_by(box,time,variable) |>
  dplyr::reframe(
    value=pomp::wquant(x=value,probs=c(0.05,0.5,0.95),weights=lik),
    name=c("lo","med","hi")
  ) |>
  ungroup() |>
  pivot_wider() -> group_traj

group_traj$pABA <- factor(group_traj$box,levels=c("05","04","03","02","01"),
                          labels=c("Control","0%","0.05%","0.5%","5%"))

Remove estimates for W for control mice

group_traj <- group_traj |> dplyr::slice(-which(group_traj$box=="05"&group_traj$variable=="W"))
group_traj |>
  filter(time<=20,variable%in%c("rbc","E","R","N","W","K")) |>
  mutate(variable=case_match(variable,
                             "K"~"Parasites",
                             "rbc"~"RBC",
                             "R"~"Reticulocytes",
                             "E"~"Erythrocytes",
                             "N"~"Indiscriminate killing",
                             "W"~"Targeted killing")) |>
  ggplot()+
  geom_line(aes(x=time,y=med,col=pABA))+
  geom_ribbon(aes(x=time,ymin=lo,ymax=hi,fill=pABA),alpha=0.2)+
  scale_y_log10()+
  scale_colour_manual(values=cbPalette)+
  scale_fill_manual(values=cbPalette)+
  xlab("Days")+ylab("Density per microlitre")+
  facet_wrap(variable~.,scales="free_y")+
  theme_bw()+
  theme(
    strip.background = element_blank()
  )
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

Plot group-level trajectories for reticulocytes with day of peak parasitaemia

Plot group-level trajectories for parasite density

## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

Plot group-level trajectories for percent reticulocytes

Create swirly data frame

## # A tibble: 6 × 29
##   box    time pABA         E      K      M      N     Qpn     Qps     Qpw    Qun
##   <chr> <int> <fct>    <dbl>  <dbl>  <dbl>  <dbl>   <dbl>   <dbl>   <dbl>  <dbl>
## 1 01        0 5%    7722711. 4.20e2 0      4.11e5 0       0       0       0.0480
## 2 01        1 5%    7939703. 2.55e3 2.55e3 4.44e5 1.52e-5 2.78e-4 1.02e-5 0.0503
## 3 01        2 5%    8165549. 1.56e4 1.56e4 4.74e5 9.34e-5 1.63e-3 1.04e-4 0.0532
## 4 01        3 5%    8166517. 9.22e4 9.27e4 4.66e5 5.29e-4 9.04e-3 1.15e-3 0.0523
## 5 01        4 5%    8085821. 5.03e5 5.18e5 4.40e5 2.45e-3 4.01e-2 1.54e-2 0.0470
## 6 01        5 5%    7717184. 2.01e6 2.29e6 4.26e5 7.74e-3 7.99e-2 1.48e-1 0.0370
## # … with 18 more variables: Qus <dbl>, R <dbl>, SM <dbl>, SN <dbl>, SW <dbl>,
## #   W <dbl>, loss <dbl>, perRetic <dbl>, rbc <dbl>, varN <dbl>, lambda_e <dbl>,
## #   lambda_n <dbl>, lambda_r <dbl>, lambda_u <dbl>, lambda_u_n_ratio <dbl>,
## #   lambda_u_w_ratio <dbl>, lambda_u_wn_ratio <dbl>, lambda_w <dbl>

Parasite density versus targeted killing

## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 7 rows containing missing values (`geom_path()`).
## Warning: Removed 2 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 3 rows containing missing values (`geom_path()`).
## Warning: Removed 3 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 14 rows containing missing values (`geom_path()`).
## Warning: Removed 2 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 21 rows containing missing values (`geom_path()`).
## Warning: Removed 5 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 21 rows containing missing values (`geom_path()`).
## Warning: Removed 5 rows containing missing values (`geom_text()`).

## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 3 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 3 rows containing missing values (`geom_path()`).
## Warning: Removed 3 rows containing missing values (`geom_text()`).

Indiscriminate killing versus reticulocyte supply

Targeted killing versus reticulocyte supply

## Warning: Removed 7 rows containing missing values (`geom_path()`).
## Warning: Removed 2 rows containing missing values (`geom_text()`).
## Warning: Removed 14 rows containing missing values (`geom_path()`).
## Warning: Removed 2 rows containing missing values (`geom_text()`).
## Warning: Removed 21 rows containing missing values (`geom_path()`).
## Warning: Removed 5 rows containing missing values (`geom_text()`).
## Warning: Removed 21 rows containing missing values (`geom_path()`).
## Warning: Removed 5 rows containing missing values (`geom_text()`).

Reticulocyte supply versus erythrocyte density

Plot N versus R with bi-directional error bars

Plot reticulocytes versus erythrocytes

Ratio of \(lambda_u\) to \(lambda_w\) over time

## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis

Ratio of \(lambda_u\) to \(lambda_w\)+\(lambda_n\) over time

## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis

Read in data of weighted individual trajectories

## Rows: 606492 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): box, mouse
## dbl (9): rep, time, E, R, W, N, K, lik, rep2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Plot individual weighted trajectories with group-level trajectories

## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 26 rows containing missing values (`geom_path()`).

## Warning: Transformation introduced infinite values in continuous y-axis
## Removed 26 rows containing missing values (`geom_path()`).